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O A MULTILAYER SAINT- VENANT SYSTEM WITH MASS EXCHANGES FOR 

SHALLOW WATER FLOWS. 
DERIVATION AND NUMERICAL VALIDATION* 

E. AuDussE^ M.O. Bristeau^, B. Perthame^'^ and J. Sainte-Marie^''' 

Abstract. The standard multilayer Saint- Venant system consists in introducing fluid layers that are 
advected by the interfacial velocities. As a consequence there is no mass exchanges between these 
layers and each layer is described by its height and its average velocity. 

Here we introduce another multilayer system with mass exchanges between the neighborhing layers 
where the unknowns are a total height of water and an average velocity per layer. We derive it from 
Navier-Stokes system with an hydrostatic pressure and prove energy and hyperbolicity properties of 
the model. We also give a kinetic interpretation leading to effective numerical schemes with positivity 
and energy properties. Numerical tests show the versatility of the approach and its ability to compute 
recirculation cases with wind forcing. 
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1. Introduction 

Due to computational issues associated with the free surface Navier-Stokes or Euler equations, the simulations 
of geophysical flows are often carried out with shallow water type models of reduced complexity. Indeed, 
for vertically averaged models such as the Saint- Venant system [7], efficient and robust numerical techniques 
(relaxation schemes [10], kinetic schemes [22],. . . ) are available and avoid to deal with moving meshes. 

Non-linear shallow water equations model the dynamics of a shallow, rotating layer of homogeneous incompress- 
ible fluid and are typically used to describe vertically averaged flows in two or three dimensional domains, in 
terms of horizontal velocity and depth variation, see Fig. 1. 

The classical Saint- Venant system [7] with viscosity and friction [14-16, 18] is well suited for the modeling 
of dam breaks or hydraulic jumps. The extended version of the Saint- Venant system proposed by Bristeau 
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Figure 1. Averaged models derived from Navier-Stokes equations. 



and Sainte-Marie [11] dropping the hydrostatic assumption is well adapted for the modeling of gravity waves 
propagation. 

Considering flows with large friction coefficients, with significant water depth or with important wind effects, 
the horizontal velocity can hardly be approximated - as in the Saint- Venant system - by a vertically constant 
velocity [24]. To drop this limitation a multilayer Saint- Venant model is often used where each layer is described 
by its own height, its own velocity and is advected by the flow (see [1,5,6] and the references therein). This 
advection property induces that there is no mass exchanges between neighborhing layers and makes a close 
relation to models for two non-miscible fluids (see [9,12,13]) for instance). In [1] the multilayer strategy was 
formally derived from the Navier-Stokes system with hydrostatic hypothesis departing from an earlier work [6] 
introducing a vertical partition of water height. 

Here, we derive another and simpler multilayer model where we prescribe the vertical discretization of the 
layers taking in to account the (unknown) total height of water. Using a Galerkin approximation in lagrangian 
formulation, we obtain a system where the only additional unknowns are the layers velocities. This leads to a 
global continuity equation and allows mass exchanges between layers. 

The objective of the paper is to present the derivation of this new multilayer model and to exhibit its main 
properties (hyperbolicity, energy equality, . . . ) . Some simulations performed with a kinetic scheme [4] are 
presented at the end of the paper. 

The paper is organized as follows. In Section 2, we flrst present, in a simplifled case, the formulation of the 
new multilayer Saint- Venant system starting from the hydrostatic Euler equations. In Section 3, we recall the 
Navier-Stokes system with a free moving boundary and its closure, and the Shallow Water system. We also 
introduce the multilayer formulation in the context of the hydrostatic assumption. In Section 4 we examine 
the main properties of the multilayer system and present a kinetic interpretation of the proposed model. This 
kinetic formulation leads to a numerical scheme detailed in Section 5 where some numerical simulations are also 
shown. 



2. A SIMPLIFIED CASE 



Before deriving the complete version of the multilayer system, we illustrate the approach in a simple situation. 
Moreover this case emphasizes the main differences with the multilayer system proposed by Audusse [1]. 
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We depart from the free surface hydrostatic Euler system 



du dw 
dx dz 
du du^ duw dp 
dt dx dz dx 
dp 

dz 



0, 
0, 



(2.1) 
(2.2) 
(2.3) 



for 

t > to, X gM., zi,{x) < z < 'i]{x,t), 
where T]{x,t) represents the free surface elevation, u — {u,w)'^ the velocity. The water height is H = i] — zi,, see 
Fig. 2. 

We add the two classical kinematic boundary conditions. At the free surface, we prescribe 



dr] dr] 



(2.4) 



where the subscript s denotes the value of the considered quantity at the free surface. At the bottom, the 
impermeability condition gives 

M&^-Wh = 0, (2.5) 
where the subscript h denotes the value of the considered quantity at the bottom. 

We consider that the flow domain is divided in the vertical direction into N layers of thickness ha with TV + 1 
interfaces 2^+1/2 (x,t), a = 0, (see Fig. 2) so that 



and 



N 



a = l 



,t) = Zb{x) + y^^hj{x,t). 



(2.6) 



(2.7) 
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Figure 2. Notations for the multilayer approach. 
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We consider the average velocities Ua, a = 1, . . . , N defined by 



1 



^1/2 



Ua{x,t) = — / u{x,z,t)dz, (2.. 



we also denote 



1 

< >a [X,t) = ■ 

and 



z„+l/2 

< >a {x,t) = I u^{x,z,t)dz, (2.9) 



Ua+1/2 u{x,Za+l/2,t), (2.10) 

the value of the velocity at the interface 2q-+i/2. 

Proposition 2.1. l^iift these notations, an integration of (2.1)-(2.3) over the layers [zq,_i/2, -2(1+1/2]; o; = 
1, leads to the following system of balance laws 

dha dKua ^ r r 1^ 

dhaUa 9 /, 2 \ J dzb ^ „n 
^ '^'^''^^ ~dx ^ -^/la"^ + "q+1/2<^q + 1/2 - Wq-1/2<-'q-1/2- (2.12) 

r/ie expression of the exchange terms Ga+\/2 given in the following. 

Proof. The proof relies on simple calculus based on the Leibniz rule. Using the incompressibility condition (2.1) 
integrated over the interval [zq„i/2, 2Q+1/2], we deduce the mass equation (2.11) where we exhibit the kinematic 
of the interface on the right hand side 

dZa+l/2 dZa+l/2 , , a AT /o 1 Q^ 

L^Q+1/2 = ^-WQ+l/2 — 2:Q+i/2,t), q; = 0, ...,iV. (2.13) 

The relation (2.13) gives the mass fiux leaving/entering the layer a. through the interface Zq,+i/2. 

Then we consider the velocity equation (2.2). We first observe that from the hydrostatic assumption (2.3) one 
can compute the pressure as a function of the water height : 

P{x,z,t) = g{r]{x,t) - z). 

Now we integrate the equation (2.2) over the interval [za-1/2, Za+1/2] and we obtain the relation 

dhaUa d 2 \ 7 (^V ^ ^ - .X 

dt 'dx ^ >aj+9ha-g^^Ua+i/2C^a+l/2-Ua-l/2(^a-l/2, (2.14) 

and with the definition of H, this is equivalent to (2.12). Then the kinematic boundary conditions (2.4) and 
(2.5) can be written 

Gi/2 = 0, Gjv+i/2 = 0. (2.15) 
These equations just express that there is no loss/supply of mass through the bottom and the free surface. 

Notice also that one can compute Ga+i/2j just adding up the equations (2.11) for j < a and using the first 
equality of (2.15) 

Gc+i/2=QiY.^, + -Q^Y.^,u,, a = \,...,N. (2.16) 

□ 
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The standard multilayer Saint- Venant system [1] is obtained by prescribing 

G„+i/2 = 0. (2.17) 

This choice is clearly natural for inmiscible fluids but is not justifled if the multilayer system is seen as a 
numerical approximation of the hydrostaic Euler equations. Indeed there is no reason to prevent the water 
exchanges between connected layers. Moreover it is exhibited in [1] that this choice may lead to the development 
of instabilities at the interfaces. 

Here we drop this assumption and we only keep the two physical kinematic boundary conditions (2.15). The 
equation (2.11) is then no nore meaningful since the quantity appears on both side of the equality. Nev- 
ertheless the sum of the equations (2.11) for all the layers is still relevant and the boundary condition (2.15) 
leads to a global continuity equation for the total water height H 

dH d 



- ^ haUa = 0, (2.18) 



dt dx 

a—l 

and each layer depth ha is then deduced from the total water height by the relation 

K = LH, (2.19) 

with la, OL = l,...,Af a given number satisfying 

AT 

> 0, ^lct = 1- (2.20) 

Thus the momentum equation (2.12) becomes 

+ ^\ha<U >a+-, 7^ ] = -ghaj— +Ua+l/2<~^a + l/2-y'a-l/2(^a-l/2- (2.21) 



dt dx \ la 2 J ' dx 

Using (2. 18), (2. 19), the expression of Ga+1/2 given by (2.16) can also be written 

J=l \ i=l / 

Finally we have to define the quantities ha < >a and Ua+1/2 appearing in (2.12). As usual in the derivation 
of such systems, we have considered ha < >q~ haU^, this will be discussed in details in paragraph 3.5. The 
velocities Ua+1/2, a = I, ■ ■ . , N — I are obtained using an upwinding 

, - / Ga+1/2 > , . 

Ua+1/2 — i „, ;f ^ ^ n 

To illustrate the formulation of the new model, we compare it with the system proposed in [1] in the simple 
case of a two-layer formulation. Neglecting the viscosity and friction, the formulation obtained by Audusse [1] 
corresponds to (2. 11), (2. 12) with (2.17), i.e. 

dhi dhiui ^ dh2 dh2U2 ^ ^^-i 
dt dx ' dt dx ' 
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dhiui 


dhTui 


dt 


1 — 

dx 


dh2U2 


1 dh2U^ 


dt 


dx 



9hiir-+u3/2[l^ + l^r-h (2.28) 



dzb 
dx 


+ "3/2 1 


[ dt 


^dHui 
dx 


dzb 
dx 


- "3/2 1 


[ dt 


^dHui 
dx 



, d(hi + /it) , dzh 

9h.^^--9h.^, (2.25) 

9h.'-^^^--9h.'^, (2.26) 

with hi + h2 = H. The preceding formulation corresponds to a superposition of two single layer Saint- Venant 
systems (see also [9, 12, 13] where a very similar model is considered in a bi-fiuid framework). 

With our approach (2. 18), (2. 21), the two-layer formulation reads 

dH dhiui dh2U2 

dhiui dhiuf g dHhi 

dt dx 2 dx 
dh2U2 dh2ul , gdHh2 . , , . . 

where hi = IH, /i2 = (1 - l)H, (2.30) 

with / g (0, 1) prescribed. The velocity at the interface, denoted M3/2, is calculated using upwinding, following 
the sign of the mass exchange between the layers. It is important to notice that, in the new formulation (2.27)- 
(2.30), we obtain directly a left hand side term written in conservative form with the topography and the mass 
exchange as source terms whereas the pressure term of (2.24)-(2.26) has to be modified [1] to get a conservative 
form. Moreover we prove in Section 4 that the system (2.27)-(2.30) is hyperbolic, which is not the case for 
system (2.24)-(2.26). 

The difference between (2.27)-(2.30) and (2.24)-(2.26) mainly comes from the physical definition of the layers. 
Audusse introduces a physical discretization where each layer has its own continuity equation. These N conti- 
nuity equations mean the layers are isolated each other, this situation corresponds to the case of N non miscible 
fluids. In the formulation (2.27)-(2.30), the discretization correponds to a finite elements approximation - of 
Po type - of the velocity u. In this case, the definition of the layers does not correspond to a physical partition 
of the flow but is related to the quality of the desired approximation over u. Thus we have only one continuity 
equation meaning the fluid can circulate form one layer to another. 

3. Derivation of the viscous multilayer shallow water system 

In this section we will apply to the Navier-Stokes equations the multilayer approach presented in the preceding 
section. 

3.1. The Navier-Stokes equations 

Let us start with the incompressible Navier-Stokes system [17] restricted to two dimensions with gravity in 
which the z axis represents the vertical direction. For simpHcity, the viscosity will be kept constant and 
isotropic throughout the paper (we refer the reader to [14] for a more general framework). Therefore we have 
the following general formulation: 



du dw 
dx dz 



0, (3.31) 

du ^ ^du ^ ^du ^ dp _ d^xx ^ d^xz (3 32) 
dt dx dz dx dx dz 

dw dw dw dp dY^^x dJ^^z ,„ 

1^ +^ir + IT = ^9 + ^ — + ~^' ^-23) 
dt ox oz oz dx dz 
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and we consider this system for 

t > to, X G R, zi,{x, t) < z < r/{x, t). 

We use the same notations as in the previous section. We now consider the bathymetry Zh can vary with respect 
to abscissa x and also with respect to time t. The chosen form of the viscosity tensor is symetric 

du ,du dw. 

ox oz ox 

dw ,du dw. 

Oz oz Ox 

with v the viscosity coefficient. 
3.2. Boundary conditions 

The system (3.31)-(3.33) is complete with boundary conditions. The outward and upward unit normals to the 
free surface Us and to the bottom nf, are given by 



dr] \ ^ / _d_ 



1 ) ^^TiM^ ' 



Let Et be the total stress tensor with 

Et = -pid 



"^xx ^xz 



At the free surface we have the kinematic boundary condition (2.4). Considering the air viscosity is negligible, 
the continuity of stresses at the free boundary imposes 

Eth,, = -p°n„ (3.34) 

where p° = p'^{x, t) is a given function corresponding to the atmospheric pressure. Relation (3.34) is equivalent 
to 

ns.ETHs = -p", and ts.Erris = 0, 

tg being orthogonal to ng. 

Since we now consider the bottom can vary with respect to time t, the kinematic boundary condition reads 

+ Ub-^ - Wb = 0, (3.35) 

where (x,t) 1-^ Zb{x,t) is a given function. Notice that the equation (3.35) reduces to a classical no-penetration 
condition (2.5) when Zb does not depend on time t. For the stresses at the bottom we consider a wall law under 
the form 

ETHfc - (nfc.ETn6)nfc = K(vb, H)wb, (3.36) 

with V(, = Ufc — (0, ^)"^ the relative velocity between the water and the bottom. If K(vb, H) is constant then 
we recover a Navier friction condition as in [16]. Introducing ki laminar and kt turbulent friction coefficients, 
we use the expression 

K{vb,H) = ki + ktH\^^b\, 
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corresponding to the boundary condition used in [18]. Another form of K(vb,H) is used in [10] and for other 
wall laws, the reader can also refer to [19]. Due to thermomechanical considerations, in the sequel we suppose 
K(vb, H) >0 and ^(vb, H) is often simply denoted by k. 

Let tf, satisfying t^.n;, = then when multiplied by t;, and rifc. Equation (3.36) leads to 

t^.S^nft ~ KVfc.tb, and v^.n^ = 0. 

3.3. The rescaled system 

The physical system is rescaled using the quantities 

• h and A, two characteristic dimensions along the z and x axis respectively, 

• Us the typical wave amplitude, ab the typical bathymetry variation, 

• C = \/gh the typical horizontal wave speed. 

Classically for the derivation of the Saint- Venant system, we introduce the small parameter 

h 

"=A- 

When considering long waves propagation, another important parameter needs be considered, namely 

and we consider for the bathymetry ^ = 0{S). Notice that e is related to a priori informations only associated 
to geometrical features whereas Qs and accordingly S deal with the state variables of the problem. 

Depending on the application, S can be considered or not as a small parameter. For finite amplitude wave 
theory and assuming Zb{x, t) = z^, one considers e ^ 1, 5 = 0(1) whereas the Boussinesq waves theory requires 

(5 < 1, e < 1 and C/^ = 0(1) 

where Ur is the Ursell number defined by = Jr, see [25]. All along this work, we consider e <^ I whereas, 
even if the parameter 6 is introduced in the rescaling, the assumption (5 ^ 1 is not considered except when 
explictly mentioned. 

As for the Saint- Venant system [16, 18], we introduce some characteristic quantities : T = A/C for the time, 
W = Qs/T = e6C for the vertical velocity, U = W/e ~ SC, for the horizontal velocity, P = for the pressure. 
This leads to the following dimensionless quantities 

X ^ Z ^ 1] 7 t 

A n, fls i 

^ p ^ u w 

Notice that the definition of the charateristic velocities implies 5 = ^ so 5 also corresponds to the Froude 
number. When 5 = 0(1) we have U ~ G and we recover the classical rescaling used for the Saint- Venant 
system. For the bathymetry z^ we write Zb{x, t) = Zb(x) + bit) and we introduce Zf, = Z^/h and b ~ bja^. This 
leads to 

dzb _ _ ^ _ ^ 9zb 

dt di or dx dx 

The different rescaling applied to the time and space derivatives of Zb means that a classical shallow water 
assumption is made concerning the space variations of the bottom profile whereas we assume the time variations 
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of Zf, lie in the framework of the Boussinesq assumption and are consistent with the rescaling applied to the 
velocity w. 

We also introduce v = ^ and we set k = ^. Notice that the definitions for the dimensionless quantities are 
consistent with the one used for the Boussinesq system [21,26]. Notice also that the rescaling used by Nwogu [20] 

2 

differs from the preceding one since Nwogu uses w = ^w. 

As in [16, 18], we suppose we are in the following asymptotic regime 

V ~ eiyQ, and k = ekq, 

with kq = Kifi + eKffii^b, H), being constant. 

This non-dimensionaHzation of the Navier-Stokes system (3.31)-(3.33) leads to 



0, 



du dw 
dx dz 

di dx 



eS' 



dw ^ duw 



dx 



o duw 
^ 


dp 


dz 


dx 


dw^\ 


dp 


^ dz 


dz 



dx 



du 



dx 



dz 



du 



dw 



dz 



dx 



dx 



du 



TTT eSi^o^ + i^o£ S— + eS— 2vq — 



d 



dx 



dz 



dw 



dz 



(3.37) 
(3.38) 

(3.39) 



where we use the divergence free condition to write velocity equations (3.38) and (3.39) in a conservative form. 
The associated boundary conditions (2.4), (3.34), (3.35) and (3.36) become 



df) 
5vo 



6u 

dw 
dz 

du 

'dS 



dfj 
dx 



= 0, 

dx \ dz 



dx 



eS 



dfj 
dx 



, dw 
dx 
du 
dx 



dx 



db dzb , ^ 
dt dx 



dx 
dzb 



du 
dz 



— [ 2eSiyo ^ 



dx 



dx 



Pb 



dx \ dz 



eSKo\ 1 + 



dzb 
dx 



Pb 



Ub 



dzb du 



dx V dz 



_2 dzb ,^ 
dx 



dh 
dt' 



e^5 — 

dx 



(3.40) 
(3.41) 

(3.42) 
(3.43) 



(3.44) 



For the sake of clarity, in the sequel we drop the symbol and we denote 4t = -§r . 



3.4. The Shallow Water system 

The derivation of multilayer approximation is somehow technical. In order to better explain the analysis we 
recall the monolayer case following the asymptotic expansion in [16]. 

In the following the two sets of equations (3.37)-(3.39) and (3.40)-(3.44) are approximated to retain only the 
high order terms. 
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Due to the applied rescaling some terms of the viscosity tensor e.g. 



dx 



dx 



are very small and could be neglected. But, as mentioned in [1, Remarks 1 and 2], the approximation of 
the viscous terms has to preserve the dissipation energy that is an essential property of the Navier-Stokes and 
averaged Navier-Stokes equations. Since we privilege this stability requirement and in order to keep a symmetric 
form of the viscosity tensor, we consider in the sequel a modified version of {3.37)-{3.39) under the form 



du dw 
dx dz 



0, 



dw „duw 



dz 

dt ' " dx ' dz 
corresponding to a viscosity tensor of the form 

du 



^yduw dp .,^9 /„ du\ d f ^ du 
dx \ dx dz \ dz 



dx 
dp 
dz 



d I ^ au\ c / ^ „ ow 
dx \ dz J dz \ dz 



du 



d 



dt 



(3.45) 
(3.46) 

(3.47) 



dx dz 



^zz = 2v 



dw 
dz 



This means the terms in e^dxW have been neglected in (3.37)-(3.39) and in (3.40)-(3.44). For details about the 
adopted form of the viscosity tensor see [11, Remark 2] and [1, Lemma 2.1]. 

In the same way, retaining only the high order terms, the boundary conditions (3.40)- (3.44) become 



dri (. dr] 
_+Su,--w, 



0, 



dw 

2e5vQ — 
dz 



5vq 



du 
dz 



2 dri du 
Ps-ed vq— — 
dx dz 



dx V dx 



dzb , dzb ^ 



6 1^0 



du 



dzb ( du 
dx V dx 



Pb 



-6p\ 



dx 



J. dw 
dx V dz 



dzb du 
dx dz 



( ^ 
I dx 



3/2 



life. 



(3.48) 

(3.49) 

(3.50) 
(3.51) 

(3.52) 



Now we will exhibit the hydrostatic and non hydrostatic parts of the pressure. An integration of (3.47) from z 
to 5-q gives 



^"^dw d{uw). 



dt 



dx 



)dz + e'^S'^iwl - w'^) + Ps -p 



\ ^ f^"^ 9 f du\ , ^ 9w „ c 
= -((577 - z) + ed J — \^o-g^ j dz- 2ediyo— + 2e5vQ — 



(3.53) 
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From the equations (3.41) and (3.42) it comes 



du 
dz 



and the boundary condition (3.41) gives 



dw 
dz 



(3.54) 



(3.55) 



The previous relation and the kinematic boundary condition (3.40) allow us to rewrite (3.53) under the form 

e'6 



w dz + S— I {uw) dz ) - e^S'^w'^ + 5p°- - p 



dt 



Classically we have 



and using relations (3.45), (3.56) and the Leibniz rule we have 

r''" d f du 



du 


^drj du 


du 


dx 


g dx dz 


s 9x 



(3.56) 



eS 



dx V '^^ dz 



, ^ c 9w ^ du ^ du 
dz - 2edi^o— = ediyQ— + edi^o-^ 
dz dx ax 



0(£3(5). 



This leads to the expression for the pressure p 

p = Ph+Pnh + 0{e^S), 
where the viscous and hydrostatic part ph is given by 

c n fc \ ^ du ^ du 

Ph = op + (djy - z) - S'J^o^ - £Ovq — 



(3.57) 



and the non-hydrostatic part Pnh is 

Pnh = e^5 



d d \ 

w dz + 5— I (uw) dz ~ e^S'^uP . 



dt 



dx 



The derivation and analysis of a classical Saint- Venant type system taking into account the non-hydrostatic 
part of the pressure has already been carried out by the authors [11]. The derivation of the multilayer system 
in this general framework is in progress. It will be presented in a forthcoming paper. 

In the sequel, we restrict to the situation Pnh ~ 0. Due to this hydrostatic assumption, we have 



and we retain for ph the expression 



p^Ph + 0{e^6), 



Ph = Sp"- + (St] - z)~ 2eSvo 



du 
dx' 



(3.58) 



(3.59) 
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Then using (3.42), (3.44) and (3.55) one obtains 



du 
dz 



du 
dz 



0{e) 



From (3. 58), (3. 59) we can write 
leading to 



(3.60) 
(3.61) 



The preceding relation inserted in (3.46) leads to 



dz^ 



and Equations (3.60) and (3.62) mean that 

u{x, z, t) = m(x, 0, t) + ©(e). 



(3.62) 



(3.63) 



i.e. we recognize the so-called "motion by slices" of the usual Saint- Venant system. If we introduce the averaged 
quantity 

u dz. 



Si] - Zb 

it is well known [11,15,16,18] that the shallow water system (3. 45), (3. 46) with an hydrostatic pressure (3. 58), (3. 59) 
is approximated in 0{e'^S) by the following Saint- Venant system written with the variables with dimension 



dH dHu 
-dt^ 
dHu 



dx 



= 0, 



5 dlP_ ^ _jjdj[f_ _ gjj^ _^ _^ UiyH—) 
2 dx dx dx dx dx 



3v 



H 



(3.64) 
(3.65) 



dt dx 
with H ^ 1] — Zb. 

3.5. The viscous multilayer Shallow Water system 

We again consider the Shallow Water system (3. 45), (3. 46) with an hydrostatic pressure (3. 58), (3. 59). Here 
another approximation is introduced concerning the velocity u, it is no more assumed constant along the 
vertical but is discretized in the z direction using piecewise constant functions, see Fig. 2. As introduced in 
Sec. 2 the interval [^b,^??] is divided into N layers of thickness ha and we use the definitions (2. 19), (2. 20). We 
write 



AT 



=(a;, z, {za}, t) = Y^ l[^„-i/2,2c.+i/2] iz)ua{x, t) 



(3.66) 



with the velocities Ua, a £ [1, . . . , A^] defined by (2.8). 

Notice that from (2.7) we have zi/2 = ^6 = C(l) and zjv+1/2 = Srj = 0{5). The difference of magnitude between 
zi/2 and 2Ar+i/2 makes the assumption (5^1 difficult to integrate in the definition of the {zq+i/2}. 

Now we try to quantify the error between u and its piecewise approximation u"*'^. First we notice that in 
absence of friction at the bottom and due to the Shallow Water assumption, the relations (3.60) become 



du 
dz 



du 
dz 



(3.67) 



TITLE WILL BE SET BY THE PUBLISHER 



13 



This means we can consider that except for the bottom layer, each layer inherits the approximation (3.67) i.e. 

— = 0{e^) for z > Z3/2, 

and therefore for all a > 1 

u{x,z,t)-ua{x,t) = 0{e'^), (3.68) 

or equivalently 

u{x,z,t)-u'"'''ix,z,{za},t) = 0{e^), for z > Z3/2. 
In the bottom layer we only have 

u{x, z, t) — ui{x, t) — 0{e), 

but as in [11,16], it can be proved that we have an approximation of the velocity through a parabolic correction 



(3.69) 



for z £ [21/2, -23/2]- Using the discretization (2. 7), (2. 8) and (3.66) we claim 
Proposition 3.1. The multilayer formulation of the Saint-Venant system defined by 

a— 1 

dhiui ^ dhiu\ ^ g dh\ ^ dp"- ^ dzf, ^ ^ 

dt dx 2li dx ^ dx ' ^ dx '^^^ 

9 / dui\ dz^du^ U2-U1 

dhaUa , dhau\ ^ g dhl _ dp'' , dzb 



dt dx 2la dx dx dx 

dzj duj 
dx dx 



+ — 4j//ia— ) -4t/ 
dx \ dx 



-ii=Q+i/2 

j=a — l/2 "ct+l + ha ha + /Iq — 1 



/orae {2,...,7V-1} 



dh^UN ^ dhpju% ^ g dh% ^ _ gh^ — 

dt dx 21m dx dx dx 



hN-^ g^N-^ UN-1/2GN-1/2 



+ # { 4-/^-^) + 4 /^"-^/^^""'^/^ - 2. r"'r"-\ (3.73) 
dx \ dx J dx dx h^ + htq-i 

with ha = laH{x,t) and Ga+1/2 given by (2.16), results from a formal asymptotic approximation in 0{£^6) 
coupled with a vertical discretization of the Navier-Stokes equations (3.37)-(3.39) with hydrostatic pressure. 

Proof. The integration of the divergence equation (3.45) on each layer has been already performed in the proof 
of Proposition 2.1. We recall that the deduced layer mass equations (2.11) are not meaningful if no hypothesis 
is made concerning the mass exchange term Ga+1/2 defined by (2.16). We thus consider a global mass equation 
(3.70) by adding them up. We can also directly integrate the divergence equation from bottom to free surface 
in order to obtain equation (3.70). 
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We now consider the horizontal velocity equation (3.46) integrated over the interval [zq_i/2, ■Za+i/2]- Using for 
each layer an approximation similar to (3. 68), (3. 69), we prove that 



v+l/2 



u^ix, z,t)dz = ul + 0{s^). 



In the context of the hydrostatic approximation, we assume that the pressure satisfies (3. 58), (3. 59). The 
treatment of the inviscid part of the pressure has already been presented in the proof of Proposition 2.1 where 
we have written for the gravitational part of the pressure 



Notice that it is also possible to write 



+ 1/2^ 

dx 



^{5ri -z)dz = ^i-hl + 



2la dx ' 



la dx 



(3.74) 



-1/2 




2q + 1/2 

dx 



N 



j=a+l 



d 



^Q-l/2 

dx 



N 



(3.75) 



The expressions (3.74) and (3.75) lead to the same property for the complete model even if the hyperbolic part 
is modified. The second formulation seems more adapted to the physical description "by layers" of the system 
but leads to complementary source terms whose discretization is subtle. We will use and analyse (3.75) in a 
forthcoming paper. In the following we use (3.74). 

The integration of the viscous part of the pressure leads to 



+ 1/2^ 

dx 



-1/2 



2v 



du 
dx 



d_ 

dx 



2vh, 



dUa 

dx 



2v 



dzj duj 
dx dx 



j=a+l/2 



]=a-l/2 



It remains to consider the viscous terms on the right hand side of (3.46). The first one is similar to the viscous 
part of the pressure term. For the second one, using finite differences along the vertical, we write 



dz 



du 



dz 



2v 



du 
dz 

+ 1/2 
Ug+l - Ug 

hg+l + ha 



2v 



du 
dz , 

^a-l/2 

Ua '^Q — 1 
ha + ha-l ' 



and relation (3.72) follows. Notice that equations (3.71) and (3.73) are concerned with the evolution of the 
discharge in the lowest and uppest layers, respectively. The difference between equations (3.71) and (3.73) and 
the general equation (3.72) comes from the particular form of the viscous effect at the bottom and at the free 
surface. 

Finally we drop the 0{e'^5) terms and recovering the variables with dimension, we obtain the system (3.70)- 
(3.73). □ 



4. Properties of the multilayer system 



In this paragraph we examine some properties of the model depicted in Proposition 3.1. We study its hyper- 
bolicity and we exhibit an energy inequality and a kinetic interpretation of the system. 
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4.1. Hyperbolicity 

For the simplicity of the discussion we mainly restrict in this subsection to the two-layer version of the multilayer 
model. 

Let us first say some words about the multilayer system (2.24)-(2.26) introduced by Audusse [1]. This non- 
miscible multilayer system was proved to be non-hyperbolic. In the general case the system exhibits complex 
eigenvalues. In the very simple case ui = U2 = u the eigenvalues of the hyperbolic part was shown to be equal to 
the classical barotropic eigenvalues of the monolayer shallow water system u + ^/gH, u — y/gH plus a baroclinic 
eigenvalue u that is concerned with the interface waves. Nevertheless the system is not hyperbolic since u is 
a double eigenvalue associated to a one-dimensional eigenspace. This lack of hyperboHcity may lead to the 
development of instabilities at the interface [1,12]. In [1] a technical trick is proposed to cure the problem. Here 
we can prove the well-posedness of the system. 

Proposition 4.1. The two-layer version of the multilayer Saint-Venant system (3.70)-(3.73) is strictly hyper- 
bolic when the total water height is strictly positive. 



Proof. The two-layer version of the multilayer system depicted in Proposition 3.1 stands - we denote u — 
with u — ui or M = U2 (see 2.23) - 



dt ^ dx ' dx 



0, 



dHui ^ dHul ^ g dH^ _ ^^dzb 
dt dx 2 dx ' dx 



dHu2 dHul 



gdH^ 



dt 



dx 



2 dx 



-gH 



dzb 
dx 



dH dHui 
dt dx 

dH dHu2 
dt dx 



The previous formulation can be written under the quasi-linear form 



with 



X = 



Six) 




Mix) = 








I 



-gH'-t + ^ ("2 - ^^i)) - ^(v, H)u, - H^ 

-aH'-§^-j^(^2-u,)-H'-£ 



AiX) = \ gH - ul 2ui - u 
gH -u\ 



(1-0 



2U2 — u 



and qi ~ Hui, i = {1,2}. 

The three eigenvalues of Af~^(X)A(X) are the roots of -D(x) = det(A — xM) ~ with 

Z?(x) = —xll^^ii2ui — u — — /(2'U2 — u — x)igH — u\ + ux) — (1 — Z)(2mi — u — x)igH — u\ + ux). 

Let us fix , /, Ul and U2 in R. Let us suppose ui < U2 with U2 = wi -I- 7^. We recall that the value of the 
interface velocity u is taken equal to ui or U2 following the direction of the exchange of mass between the two 
layers. 
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Let US first suppose that u = ui. Then we obviously have 

D{ui) ^ -2gHlj^ < 0, D(-oo) = +oo, D(+oo) = -oo, 
and some computations lead to 

+ 2^72)) > 0, 

since D{u2) = (1 - 20,9^7^ + ^7^ > if ^ < 1/2 and D{ui+2l-i^) = 2^(1 - /)(4Z - 1)7^ > if / > 1/2. It follows 
that D{x) has three real and simple eigenvalues. 

Let us now suppose that u = U2. Then we have 

D{u2)^2gH(l-l)j'^ >0, D(-cx)) = +00, D{+oo) = -00, 
and some computations lead to 

D{min{ui =U2- 7^ 1*2 - 2(1 - 07^)) < 0, 

since D{ui) = {1 -2l)gH-f^ - {1 - l)-/^ < if ^ > 1/2 and i:»(w2 - 2(1 - ^7^) = 2?(1 - 0(4/ - 3)7*^ < if / < 1/2. 
Here also D{x) has three real and simple eigenvalues. 

The case U2 < is similar and we can conclude that the two-layer version of the multilayer system depicted in 
Proposition 3.1 is strictly hyperbolic. Notice that when ui = U2 = u, we find the same baroclinic and barotropic 
eigenvalues u, u + y/gH, u — \/gH as for the nonmiscible multilayer system [1], but they are all simple eigenvalues 
in this case since we consider a system with only three equations. □ 

In the case of N layers the matrices A{X) and M {X) can be written 



A 



N+l 



( h ... 

gH - u\ 2ui vi^2 

■ V2A 



\ gH - u% VN,i 
with Vij = Ui_i/2 * Ij/k and Vij = Ui+1/2 * h/li, 

( 1 

VI 1 

7V+1 = : 



In \ 



VN-l.N 

2uN J 



M 



\ 





1 / 



V 

with iH = u,_i/2 * hl^i + ""4+1/2 * hl^^■ 

We have perfomed various numerical evaluations of the eigenelements of the matrix M^^j^^An+i with numerous 
choices of the parameters H, Ua, Ua+1/2 and la- All these tests have always shown that the matrix is diag- 
onalizable on M. In the simple case where all the layers have the same velocity u, the barotropic eigenvalues 
u + \/gH and u — \fgH are simple and the baroclinic eigenvalue u has a multiplicity of iV — 1 but the matrix 
remains diagonalizable on R and the problem is still well-posed. 
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4.2. Energy equality 

The classical Saint- Venant system (3.64)-(3.65) admits an energy equality [1,11] under the form 



dt dx 



di 



du^2 k{-v,H) 



dt 



dx 



3i/ 



dt 



with Esv = ^ + aHiv+^b) _^ jj^a jjgj.g Yia^e the following result 

Proposition 4.2. For i/ie multilayer Saint-Venant system (3.70)-(3.73), smooth solutions satisfy the energy 
equality 



d_ 

di 



at 



dt 



with ET,"^ 



Proof. The proof rehes on classical computations. Starting from (3.46) with u = u™"^ , p = ph multiplying it 
with u™'^ and integrating over [za-1/2, ■Za+i/2] with 1 < a < we obtain 



dt ^"^^ ax 



-1/2 



^^+-„-l/2^^--«-l/2 



-l'OUa-1/2 



5z 



/2, 



ZC.-1/2 
Q^rnc 



dz 



2c + l/2 



5x 



1^0 



5Za+l/2 

dx 



Q+1/2 



- "a+1/2 



dt 



1+1/2 



9x 



2a + l/2 



dufy 2 dp^ dzfj 



where we have considered for z E [Za-1/21 ^a+1/2] 



dt 



du 1 , 

9^ = 7^(""+i/2~""-i/^^ 



(4.78) 



An analoguous calculation is valid for a = 1 and a = TV. A sum from a = 1 to a = TV of the equalities (4.78) 
with the boundary conditions (3.48)-(3.52) completes the proof. □ 



4.3. Kinetic interpretation 

For the simulation of a multilayer system several strategies are possible. Pares et al. [13] consider the full 
system and build a specific solver for the two-layer case. Following the discrete multilayer scheme proposed by 
Audusse [1] we prefer to exhibit a kinetic formulation of the system obtained in Proposition 3.1. Indeed kinetic 
schemes might be one of the best compromise between accuracy, stability and efficiency for the resolution of 
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Saint- Venant type equations, see [4,22]. We refer to the next section for the presentation of the numerical 
scheme. Here we focus on the kinetic interpretation of the system. 

The kinetic approach consists in using a description of the microscopic behavior of the system. In this method, 
fictitious particles are introduced and the equations are considered at the microscopic scale, where no discon- 
tinuities occur. The process to obtain the kinetic interpretation of the multilayer model is similar to the one 
used in [4] for the monolayer shallow water system. For a given layer a, a distribution function Ma{x,t,^) 
of fictitious particles with microscopic velocity ^ is introduced to obtain a Hnear microscopic kinetic equation 
equivalent to the macroscopic model presented in proposition 3.1. 

Let us introduce a real function x defined on M, compactly supported and which have the following properties 

Xi-w) = x{w) > 
ImXiw) dw = J^w^xiw) dw ^ 1. 

Now let us construct a density of particles Ma{x,t,^) defined by a Gibbs equilibrium: the microscopic density 
of particles present at time t in the layers a, in the vicinity Aa; of the abscissa x and with velocity ^ given by 

M„(x, t, - L^^^X (^i^^^^) , a = l,...,N, (4.80) 

with 

2 gH 

Likewise, we define Na+i/2[x,t,S) by 

Na+i/2{x,t,i) =Go,+i/2{x,t) 5{S,~Ua+i/2{x,t)) , a = 0, ...,iV, (4.81) 

where i5 denotes the Dirac distribution. The quantities Ga+i/2, Q < a < N represent the mass exchanges 
between layers a and a + 1, they are defined in (2.16) and satisfy the conditions (2.15), so N112 and 
also satisfy 

^1/2(2;, t, i) = iVw+i/2(x, t, - 0. (4.82) 
We also introduce the densities Ma{x,t,S) that will be used for the energy equations , they are defined by 

( , 5^^(M)Ma^0 f£,-Uaix,t) 

IVla(X, 1,1;) — X 

Notice that the introduction of this second family of densities is not needed when we consider the two dimensional 
shallow water system. Here they take into account some kind of transversal effect at the kinetic level that is 
implicitely included into the macroscopic one dimensional shallow water system. We refer the reader to [4, 23] 
for more details. 

With the previous definitions, dropping the viscous, and friction terms, we write a kinetic representation of the 
multilayer Saint- Venant system described in proposition 3.1 and we have the following proposition: 

Proposition 4.3. The functions (if, u""^) are strong solutions of the multilayer Saint-Venant system (3.70)- 
(3.73) if and only if the set of equilibria {Ma{x,t,£^)}^^i is solution of the kinetic equations 



= Qc.{x,t,0, (4.83) 
a = l,...,7V, 
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with {^^0+1/2(2;, t, 01^=0 satisfying (4-81), (4-82). The set of equations (4-83) can also be written under the 
form 

= E + - ^ (P'' + ^ - Q.) , a = 1, . . . , iV. (4.84) 

The quantities Qa{x,t,£^) are "collision terms" equals to zero at the macroscopic level i.e. which satisfy for a.e. 
values of (x, t) 



QadS, = 0, / = 0. 

The solution of (4. 83), (4. 84) is an entropy solution if additionally 



^+^^ = Q.ixXO. a = l,...,N, (4.85) 

with ^ 

y Oa + Qo^ d( < 0. 

Proof. As previously we denote X = {H, qi, . . . , qN)'^ the vector of unknowns with qa = laHua. We introduce 
M = [Ml,. . . ,MnY and an (iV + 1) x TV matrix /C(^) defined by /Cij = 1, /Ci+ij = 5,^ ^ with 5ij the 
Kronecker symbol. 

Using the definition (4.80) and the properties of the function x, we have 

lo,H{x,t)= [ M^{x,t,Od^, (4.86) 



and 



X{x,t) = J ICi£_) M{x,t,^)d^. (4.87) 



The proof is obtained by a simple integration in ^ of the set of equations (4.83) against the matrix /C(^). First, 
an integration in ^ of (4.83) gives the continuity equation (2.11) i.e. 

dlaH dlaHUa 

and by summation we have (3.70). Actually from the definition (4.81) of Na+1/2 we have 

Na + l/2{x,t,Cldi = Ga + l/2{x,t), 



and 



/ £,Na+i/2{x,t,£,)d^ ^ Ua+l/2Ga+l/2 
JR 



Likewise for the energy balance of the layer a we proceed an integration in ^ of (4.83) against ^^/2. Since we 
have 



+ M^j d( = -fui + -^h^H, (4.88) 
C ( |-Af„ + M^^ d( = ^ul + gh^Huo,, (4.89) 
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and for the source term 



we obtain the equahty 



d f ha 2 9 1 / \ 1 c 



1+1/2, 



-~-Q^{{gZb+P°')haUa) - {gZb+p°-) 

- (ffZfc + p") G^+i/a + {gzb + Gq_i/2, 



dx 



dp" dzb 

/2 —^.+1/2 ah^-Q^ 



(4.90) 



(4.91) 



The previous relation corresponds to (4.78) where the viscous and friction terms are neglected. The sum of the 
equations (4.91) gives the energy equality for the global system and that completes the proof. □ 



The formulation (4.83) reduces the nonlinear multilayer Saint- Venant system to a linear transport system on 
nonhnear quantities {Ma}a=i^ {-^Q+i/2}a=o for which it is easier to find a simple numerical scheme with good 
theoretical properties. In the case of a single layer, for a detailed proof of the kinetic interpretation refer to [4] 
and for the treatment of the source term at this microscopic level see [23]. Notice that the choice of the function 
X remains quite open at this stage since several functions satisfy the requested properties. Following this choice 
the deduced kinetic scheme will have different properties. 



5. Numerical results 

In the applications discussed here, we assume p° = and we neglect the horizontal viscosity. Then the + 1 
equations of the multilayer system (3.70)-(3.73) can be written with the general form 

a—l 

d{laHUa) d 2 : 9 , tt2^ , rj^^b , ^ ^ 
dt 9x 2 ^ + Ua + l/2'~^a+l/2 — "0-1/2^^0-1/2 



with 



k{u, H) if a = 1 
if a^l 



— Ka{u,H)ua, a = l,...,N, (5.93) 



if a = 

ly if a = l,...,N -1 

if a = N 



The previous system is of the form: 

dX dF{X) 



, - Sb{X) + SeiX) + Sv{X) (5.94) 
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with F{X) the flux of the hyperbolic part, Sb{X) the topography source term, Se{X) the mass transfer source 
term and Sv(X) the viscous and friction terms. 

To approximate the solution of the multilayer Saint- Venant system, we use a finite volume framework. We 
assume that the computational domain is discretised by / nodes Xi. We denote Ci the cell of length Axi = 
Xi+i/2 — 2^1-1/2 with Xi+if2 = {xi + Xi+i)/2. For the time discretization, we denote = J2k<n^^'' where 
the time steps At*^ will be precised later though a CFL condition. We denote X" ~ (iJ", .. .,q^^) the 
approximate solution at time on the cell Ci with ^ = laH^u^ ^. 



5.1. Time discretization 

For the time discretization, we apply time splitting to the equation (5.94) and we write 

+ - Sb{X-) + SeiX-), (5.95) 

vn-\-l vn-\-l 

— 5w(X",X"+i) =0. (5.96) 

Classically we first compute the hyperbolic part (5.95) of the multilayer system by an explicit scheme. This first 
computation includes the topographic source term in order to preserve relevant equilibria [2] and also defines 
the mass transfer terms. Concerning the viscous and friction terms (5.96) that are dissipative, we prefer a 
semi-impHcit scheme for reasons of stability. 

5.2. Numerical scheme : explicit part 

To perform the explicit step we deduce a finite volume kinetic scheme from the previous kinetic interpretation 
of the multilayer system. Notice that even if the system is hyperbolic, the eigenvalues are unknown. Thus any 
solver requiring the knowledge of the eigenvalues while but the kinetic scheme is easily extended [5] . 

Starting from a piecewise constant approximation of the initial data, the general form of a finite volume method 
is 



j^n+l _ xj 



T^n 17111 



= AeSb'l + AeSel, (5.97) 



where ct" = At^ / Axi is the ratio between space and time steps and the numerical fiux F"^^^^ is an approximation 
of the exact fiux estimated at point a:i_(.i/2- 

The topographic source term Sbf is not deduced from the kinetic interpretation (see [23]) but computed by 
hydrostatic reconstruction, see prop. 5.1. As in [4,6] the kinetic interpretation (4.83) is used to precise the 
expression of the fiuxes F^^^^ in (5.97). First, by analogy with (4.80) we define the discrete densities of 
particles Af^ j by 

MZAO'L^xC-^). wither../^ 



Then the equation (4.83) without the atmospheric pressure and topographic terms is discretised for each a by 
applying a simple upwind scheme for the advection term 

(0 = MiM) - ^< - A^:.-i/2(e)) + At- (c+'/t.io - K'^l'/U^)) , (5.98) 

where 
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and the terms Nl''tlin . will be defined in the following 



Q + l/2/j 

We define the vectors f^+\0 = (/"^(O, • ■ • , fN^iOV, M^O = (M^O, M^^^{i)Y . Each new density 
function /^^^ is not an equilibrium but thanks to the property of the right hand side of (4.83), by analogy with 
(4. 86), (4. 87) we can recover the macroscopic quantities at time We write 



(5.99) 



and by a simple integration in ^ of (5.98) against /C(^), we can precise the macroscopic formula (5.97) (without 
the topographic term) 



If we denote 
we define 



e/C(0 M^O d^- 



More precisely the expression of F^{Xi) can be written 



( PH{X^) \ 



with 



i^+(XO = ^F+(X,) = ^Z«i? / 

a=l a=l ■'^ 

F+^{Xi)^l^H [ 

Jw>- 



{Ua,i + Wq)x(w) dw, 



We denote also 



+ wct)'^xiw) dw. 



Fh^,i+i/2 - J^/.„,,-i/2 - F+^ (X,) + F-^ (X,+i) - (F+ (X,_i) + F-^ [Xi)) 



(5.100) 



(5.101) 



(5.102) 



(5.103) 



This kinetic method is interesting because it gives a very simple and natural way to propose a numerical 
fiux through the kinetic interpretation. If we can perform analytically the integration in (5.102), i.e. if the 
probability function x defined in (4.79) is chosen to be simple enough, it is also numerically powerfull because 
the kinetic level disappears and the scheme is written directly as a macroscopic scheme for which only very 
simple computations are needed. In this paper we have used 




Let us now precise the terms N^^'^^^'^ . and so the exchange terms 5*6" defined by 



Se7 = / m) {Ktl/UO - K-ljlm) d^. (5.104) 
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From the conditions (4.82) we prescribe 

KnfiO = NlX%^^{0 = 0. (5.105) 

So we recover •S'e^^ = and the equation (5.100) defines . By summation of (5.98) we have 

a 

^t"Ktl/lM) = E - ^^fM) + ^< {K^+iMO - M-^-i/2iO)) , a = l,...,N-l, (5.106) 



and we define 
so we can write 



G:X% = / Ktl/U^)d^' a = 0,...,N, (5.107) 



a 

^t"G:l\%^^ = E N^"^' - + <(^r,.+i/2 - Fhj,-i/2)\ , a^l,...,N. (5.108) 



Then using the discrete mass conservation equation giving H^'^^, the terms G^'^\^^'^ - can be written under an 
explicit form (see (2.22)) i.e. depending only of X" 



N 



ao^.g::;;!. = e k^,^ - e ' (5-109) 
j=i \ p=i ) 

we have to notice that this definition is compatible with the free surface condition of (5.105). 
We define 

- ^z'lk ^ - <+i/2.) > (5-110) 

with, according to (2.23) 

f u" , 1 , if G"iy,^ - > 0, 

t ^a.i II '-^Q + l/2,i ^ ^• 

Then the exchange term 5'ef in (5.104) is completely defined. 

We have denoted the approximations in time of N^j^-^j^ and Gq,_(.i/2 with an upperscript n + 1/2 because we 
have to define H"^^^ at the macroscopic level to obtain the microscopic approximation of Na+ijiA which is used 
for the computation of the momentum laH'^^^u^^ . 

The source term S'&f = (56^ j, 5'6" j, . . . , ,) is an approximation of the topographic source terms. For 
stability purpose, see [2] we use the following discretization 

56^,, = 0, 56S, = /a(|(i?r+i/2-)'-f(ffr-i/2+)') (5.111) 

with 

■26,i+i/2 = max{zfc^i,Zfc^i+i}, 

= + ^h.^ - ^6,,;+l/2, (5.112) 

And we have the following proposition 
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Proposition 5.1. The discretization of the source terms given by (5. Ill), (5. 112) preserves the steady states 

{<,. = 0}^=i, Hl' + Zb,,=CsteR, Vi, yn. 

given by a "lake at rest". 

Proof. For the proof of this proposition, the readers can refer to [2]. □ 
The scheme explained in this paragraph allows to calculate given by (5.95) and (5.97). 

5.3. Numerical scheme : implicit part 

Now we aim to calculate X"'^^ from (5.96). Neglecting the horizontal viscosity, the vertical viscosity source 
term can be interpreted as a friction term between one layer and the two adjacent ones. As usual we treat this 
friction term implicitly. This leads to solve a linear system. 

The implicit step does not affect the discrete water height therefore 

and the computation of the new velocities {u2'\^}a=i leads to solve a tridiagonal N x N linear system that 
reads 

rpn,n+ljjn+l _ ^" + 1 

With C/f+i = «+\ . . . , u-+ir, q^' = m\ rk-V and 



Tr^\a,a) = l^Hr' + ^f^^ + ^^), for a G {2, . . . , TV}, 

Tr^\a,a + 1) = -^f^^), for«G{l,...,iV-l}, 

rr'-'ia-l^a) = -^{j^^Tlh)' iorae{2,...,N}. 
For the friction at the bottom, several models can be used among which are Navier, Chezy and Strickler laws. 
5.4. Stability of the scheme 

We now establish the stability property of the kinetic scheme. Classically for the Saint- Venant system, a CFL 
condition ensures the water height is non negative. This CFL condition means that the quantity of water leaving 
a given cell during a time step At" is less than the actual water in the cell. 

For the multilayer Saint- Venant system we have the same kind of requirement concerning the time step At" . 
But due to the vertical discretization, the water can leave the cell Ci of the layer a either by the boundaries 
Xi±i/2 or by the interfaces Za±i/2, see Fig. 3. This makes the CFL condition more restrictive and we have the 
following proposition 

Proposition 5.2. Assume that the function x has a compact support of length 2wm then under the CFL 
condition 

At"< min min l^H^^x, ^ ^^^^ 



l<a<N iel /, , \ . 



^a-l/2,i 
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Figure 3. Interpretation of the CFL condition for the classical Saint- Venant system (up) and 
for the multilayer system (down). 



the kinetic scheme (5.97), (5.111) and (5.102) keeps the water height positive i.e. H" > if it is true initially. 

dx 



Notice that this condition does not depend on 



Proof. The proof has been adapted from those given in [3,23]. To prove the stability property of the scheme, 
we come back to the kinetic interpretation and we proceed by induction. We assume that if" > 0, Vi and we 
prove that > 0, Vi. 

From the definition of the functions Ma in (4.80) and the positivity of the function x, we deduce 

M^_i>0, Vi, for a = 1,..., TV. 

We now introduce the quantities 

K]+ = max(0,C), K]- =max(0,-O, 
and so we can write the upwind microscopic scheme (5.98) 



= (1 - am) Ml, + amUM2,.-i + '^rKl-M^.^+i 



At" 



N 



n+l/2 
a+l/2,i 



N 



ra+l/2 
.+1/2,1 



N 



n+l/2 
'-1/2,1 



N' 



n+l/2 
Q-l/2,i 



The quantity 



N 



n+l/2 
a+l/2,i 



N 



n+l/2 
a-l/2S 



(5.114) 



26 



TITLE WILL BE SET BY THE PUBLISHER 



represents, at the microscopic level, the water leaving the cell Ci of the layer a during At". A sufficient condition 
to obtain the stability property, i.e. 



lc.H^^^= [ f^+^d^>Q, Vi, fora = l 

JR 



, . . . , ^ . , 



is then 

and this requirement is satisfied when 



N' 



-1/2 
Q+l/2,i 



N' 



n+l/2 
Q-l/2,i 



G 



1+1/2 
a+l/2,i 



G 



1+1/2 
Q-l/2,i 



We recall that we have obtained in (5.109) an explicit form of G^^Y^"^ ^. If At" satisfies (5.113), then the 
condition (5.115) is satisfied and that completes the proof. □ 



(5.115) 



(5.116) 



5.5. Second order scheme 



The second-order accuracy in time is usually recovered by the Heun method [8] that is a slight modification of 
the second order Runge-Kutta method. The advantage of the Heun scheme is that it preserves the invariant 
domains without any additional limitation on the CFL. 

We also apply a formally second order scheme in space by a limited reconstruction of the variables. An advantage 
of the new multilayer approach with only one continuity equation is that the water height can be reconstructed 
while preserving the mass conservation without difficulty. 



5.6. Numerical simulations 
5.6.1. Transcritical flow over a bump 

We first consider an academic test case that is very commonly used for the vaHdation of classical one-layer 
shallow water solvers. Here we add some friction at the bottom in order to compare solutions of one-layer and 
multilayer shallow water systems with the solution of hydrostatic incompressible Navier-Stokes equations. We 
impose an infiow (left boundary) of 1.0 m^.s~^ and the water height at the exit (right boundary) is prescribed 
to be equal to 0.6 m. The Strickler friction coefficient at the bottom is 30 to^/'^.s"^ and the kinematic viscosity 
is 0.01 m^.s~^. The data are chosen such that the ffow is supposed to reach a stationary regime that presents 
some transitions between sub- and supercritical parts and an hydraulic jump. Notice that an analytical solution 
exists for this test in the case of a single layer [4, 23]. 

The simulation results are depicted in Fig. 4, 5 and 6. The presented results correspond to a time instant tf 
where the permanent regime is achieved. Notice that we present some results related to the vertical velocity in 
Fig. 5. Since we consider a shallow water type system we do not need this vertical velocity for the computation. 
But it is possible to recover it for postprocessing purpose : departing from the computed horizontal velocity we 
use the divergence free condition (3.31) and the non penetration condition at the bottom (2.5) to evaluate an 
approximation of the vertical velocity. Notice also that the actual computations are purely one dimensional. 
Hence Fig. 4 and 5 present velocity results on a postprocessing mesh that is constructed departing from the Id 
mesh by the use of the computed layer water heights. 

The results depicted in Fig. 4 and 5 are consistent with computations performed using the hydrostatic Navier- 
Stokes equations [6] and also using the former multilayer Saint- Venant system [1] . The results depicted in Fig. 6 
exhibit that the presented solver is quite robust since it is able to compute transcritical solutions and shock 
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waves even when a large number of layers are considered. Notice also that the hydraulic jump appears to be 
overestimated by the one-layer computation when compared with other results - see Fig. 6. 




Figure 4. Horizontal velocities {ua{x,tf)}^^i with iV = 15 layers. 



5.6.2. Wind effects 

We claim in the introduction that the great interest of the new multilayer formulation that we proposed here 
is to allow mass exchanges between layers. This effect is exhibited in the numerical test that we present now. 
We consider a lake with a non trivial bottom and vertical shores. We impose a constant wind stress (from left 
to right) at the free surface. The flow is then supposed to reach a stationnary state that includes some water 
recirculations in the lake. Notice that this kind of stationnary flows is clearly impossible to compute with the 
classical one-layer shallow water system since the velocity is imposed to be constant along the vertical. They 
are also out of the domain of application of the former multilayer shallow water system that was introduced by 
Audusse [6] since they clearly involve large mass transfers (at least near the shores) between the layers. 

As for the previous case we use a reconstruction strategy in order to estimate a vertical velocity fleld and we 
present the results on a postprocessing 2D mesh that is presented in Fig. 7. In Fig. 8 we present the two 
dimensional velocity vectors on this 2D mesh. The results exhibit a global recirculation that is combined with 
two local recirculations that are induced by the topography of the lake. The qualitative aspect of the solution 
is consistent with the previsions. 
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Figure 5. Vertictal velocity {wa{x,tf)}^^i with TV = 15 layers. 
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Figure 6. Shape of the free surface for simulations carried out with different number of layers. 
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Figure 7. The geometrical model with the horizontal mesh and the vertical discretization by layers. 




Figure 8. A wind blow from the left part of the domain to the right part. The arrows represent 
the velocity field in the lake. 



6. Conclusion 

In this paper, the authors have described an exchanging mass multilayer Saint- Venant system. The derivation 
of the model, the study of its main properties and a numerical scheme for its discretization are given. Some 
simulations are also presented. Notice that the model and the results presented here in 2D (x, z) are also 
available in 3D {x,y,z). 

Because of its accuracy and simpHcity, the kinetic scheme seems well adapted for the simulations of such a 
model. Moreover since the eigenvalues of the hyperbolic system are not explictly known, a lot of finite volume 
schemes fails in this situation. 

As depicted in Fig. 7, the vertical discretization proposed for water height leads to a regular mesh. A strategy 
of "mesh refinement" based on a inhomogenous number of layers have to be added. 

The presented system can be enriched in several ways. First, the hydrostatic assumption concerning the pressure 
terms can be relaxed leading to the models presented in [11]. Then we can also consider a passive pollutant in 
the fiow. This implies to add a conservation equation for the pollutant concentration. Finally, we can consider 
the density of the fiuid varies with the concentration of pollutant. These three improvements have been added 
to their model by the authors and will be presented in forthcoming papers. 
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